Identifying locations susceptible to micro-anatomical reentry using a spatial network representation of atrial fibre maps

Micro-anatomical reentry has been identified as a potential driver of atrial fibrillation (AF). In this paper, we introduce a novel computational method which aims to identify which atrial regions are most susceptible to micro-reentry. The approach, which considers the structural basis for micro-reentry only, is based on the premise that the accumulation of electrically insulating interstitial fibrosis can be modelled by simulating percolation-like phenomena on spatial networks. Our results suggest that at high coupling, where micro-reentry is rare, the micro-reentrant substrate is highly clustered in areas where the atrial walls are thin and have convex wall morphology, likely facilitating localised treatment via ablation. However, as transverse connections between fibres are removed, mimicking the accumulation of interstitial fibrosis, the substrate becomes less spatially clustered, and the bias to forming in thin, convex regions of the atria is reduced, possibly restricting the efficacy of localised ablation. Comparing our algorithm on image-based models with and without atrial fibre structure, we find that strong longitudinal fibre coupling can suppress the micro-reentrant substrate, whereas regions with disordered fibre orientations have an enhanced risk of micro-reentry. With further development, these methods may be useful for modelling the temporal development of the fibrotic substrate on an individualised basis.


Introduction
Atrial fibrillation (AF) is the most common cardiac arrhythmia with significant impacts on both morbidity and mortality [1]. Despite extensive research, there are significant disagreements within the cardiac electrophysiology community as to the mechanisms underlying AF [2,3].
Most studies into the mechanistic origin of AF focus on the maintenance of AF, typically arguing for either organised (mother waves or stable rotors) or disorganised mechanisms a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 (multiple reentrant wavelets). However, if both organised and disorganised mechanisms of AF coexist on a continuous spectrum of electromechanical organisation, as suggested recently [4], it is unlikely that any single treatment strategy will be successful across the full spectrum of AF mechanisms. It is for this reason that personalised, patient-specific approaches to AF treatment have become a key research focus in recent years [5,6].
Here we computationally investigate one mechanism of AF initiation and maintenance: the formation of micro-anatomical reentrant circuits, which are continuously activated electrical circuits anchored to the fibre structure of the atria [7]. This mechanism is controversial, with supporting clinical evidence still required, but it has some "potential to unify . . . previous discrepant observations" in the AF literature [8]. Importantly, the size of experimentally identified micro-reentrant circuits is often at, or below, the spatial resolution which can be resolved with conventional multi-electrode mapping [7,9]. Similarly, the interstitial fibrosis which insulates these circuits is often not easy to detect using conventional LGE-MRI [4], and extracting precise fibrosis densities is challenging given the variability in signal thresholding choices between scans and patients [10]. Hence, computational approaches allow for a degree of hypothesis testing which avoids some of these challenges.
In this proof of concept, we introduce a novel method to assess the feasibility of patient-specific predictions for the distribution of micro-anatomical reentry across the atria. In particular, the method is inspired by the idea that it may be possible to predict the emergence of microreentry before the atria have accumulated sufficient interstitial fibrosis, in contrast to most other approaches which take a static view of the AF substrate [5,11].
Starting from image-based models of the atria, we combine data regarding the atrial geometry and the underlying myocardial fibre structure to form a spatial network. By progressively removing connections in the atrial structure, we assess where in the network a path exists which is sufficiently long to harbour micro-reentry, an approach related to the study of percolation in network science [12]. In the current work, our primary goal is to understand the utility of the method and demonstrate its future potential. With this in mind, our research focuses on how the regions predicted by our method of being susceptible to micro-reentry depend on fibre structure and atrial geometry.

Aims
Experimental evidence from explanted human hearts has suggested that "AF may be driven by microanatomic reentrant AF drivers anchored to fibrotically insulated tracks within the complex atrial wall" [9]. To ensure that a given fibre tract can sustain micro-reentry, the reentrant pathway through a fibrotically insulated fibre tract must be at least one refractory wavelength long.
Our approach is illustrated schematically in Fig 1. Denoting the minimum refractory wavelength during fibrillation as τ, our aim is to assess where in the atrial structure fibrotically insulated reentrant pathways exist of length ℓ, such that ℓ � τ. In a healthy heart with minimal interstitial fibrosis, atrial myofibres are well coupled in both the longitudinal and transverse directions. Hence, the longest possible reentrant path is significantly shorter than one refractory wavelength, ℓ � τ. As the density of electrically insulating interstitial fibrosis increases, the longest reentrant pathway increases until a threshold density is reached at which ℓ � τ and the fibrotically insulated fibre tract can sustain continuous reentry during AF. Importantly however, this threshold density may vary across different regions of the atria; some regions may be susceptible to micro-reentry at a low fibrosis density, whereas other regions may require a significantly higher fibrosis density before micro-reentry can be induced. If the density of fibrosis is too high in a region, micro-reentry may be prevented due to the absence of any closed loops (compact fibrosis).
Our hypothesis is that the process of electrically insulating fibre tracts with interstitial fibrosis can be modelled using spatial networks. A spatial network is a graph consisting of nodes representing entities in space, and edges representing spatial connections between those entities [13]. In our case, a node represents one or more atrial cardiomyocytes in a specific atrial region, and an edge between two nodes indicates whether those groups of myocytes are electrically coupled, that is whether an activation wavefront can pass directly from the group of myocytes represented by one node, to the group of myocytes at the other node.
Starting from a network in which the local density of edges is high, our assertion is that the accumulation of insulating interstitial fibrosis in the atria is structurally equivalent to the progressive removal of edges in our spatial network. If the density of fibrosis is low, then the density of edges in the spatial network is high, and similarly, if the density of fibrosis is high, the In the healthy atrial myocardium, cardiomyocytes (pink cells) are well coupled along and across the principle fibre directions, such that the longest reentrant loop is too short to sustain micro-anatomical reentry (green arrows). (B) Over time, the atria accumulate interstitial fibrosis (yellow cells) forming short segments of electrically isolated fibre. (C) The length of these segments grows as more interstitial fibrosis accumulates until the isolated segments are sufficiently long to harbour micro-reentry (blue arrows). We hypothesise that the accumulation of interstitial fibrosis can be modelled as a spatial network. Right column: Nodes (yellow points) representing a group of cardiomyocytes are connected to their neighbours along fibres (red links) and across fibres (blue links). If the density of interstitial fibrosis is low, there are many transverse (blue) links. We model the increase in the density of interstitial fibrosis as the progressive removal of transverse links, and equate density of edges in the spatial network is low. Then, by identifying loops in the spatial network of length ℓ � τ, we can identify which atrial regions may act as a substrate for micro-anatomical reentry at a given fibrosis density.

Materials and methods
In this paper, we generate spatial networks from three atrial datasets (see below), each consisting of a voxel mesh with local fibre orientation vectors. For each dataset, we generate spatial networks according to two different approaches: one which preserves atrial fibre orientations in the spatial network (the fibre model), and a null model approach which ignores the fibre orientation data but retains other structural information such as wall thickness (the fibre-less null model). These two approaches are shown schematically in Fig 2. By utilising both approaches, we are able to identify how the underlying atrial fibre structure affects the observed microreentrant substrate.

Atrial datasets
The spatial network algorithm is tested on three atrial datasets, one human dataset with MRI derived geometry, but synthetic fibre orientations, and two sheep atrial datasets where the fibre structure has been inferred from high resolution serial surface imaging. Each dataset provides a 3d grid of vectors indicating the local fibre orientation at each point in the atria. The sheep datasets are high resolution and inferred directly from the anatomy of each dataset. Summary statistics for the three image-based atrial fibre maps are given in Table 1. 3.1.1 Synthetic atrial fibre map. The human dataset is detailed in [15] based on techniques developed in [14]. The atrial geometry used provides no information as to the atrial fibre structure. Therefore, a synthetic fibre map is constructed using a semi-automated rulesbased approach which is described in the supplementary material (SM) section 1.1.1.
3.1.2 Anatomically derived fibre maps. We study two anatomically derived sheep atrial datasets: one healthy sheep acquired in [16] and a sheep with pacing-induced heart failure [17]. These datasets are individual-specific, retaining significant local heterogeneity in the fibre maps, and are provided at high resolution, in contrast to the synthetic human fibre map.
Information regarding the data acquisition process is provided in SM section 1.1.2. In this proof of concept, our aim is not to assess the differences in micro-anatomical reentry between humans and sheep, or between healthy and heart failure (HF) sheep, but to understand how different sources of fibre orientation data affect our results.

Constructing spatial networks with fibre structure
The schematic shown in Fig 3 outlines the process by which each image-based model is converted into a spatial network in which the underlying fibre structure is retained.
3.2.1 Fibre tractography. Fibre tracts are generated across the atrial geometry by applying a modified version of the Evenly Spaced Streamlines (ESS) algorithm, see [18]: (1) Local fibre orientation data for a given atrial dataset is coarse grained to a specific voxel resolution, . The ESS algorithm was modified to maximise fibre length by generating longer fibres first. Importantly, unlike other common tractography methods, this modified ESS approach ensures that the density of fibres is approximately uniform across the atrial structure. Technical details for the tractography process are provided in SM section 1.2.
3.2.2 Generation of spatial network structure. The fibres generated from tractography are not connected to each other, Fig 3(C). To connect fibres in a single spatial network, nodes are placed at even intervals (steps of 1 voxel length) along fibres starting at a random seed point. These nodes connect to neighbouring nodes in the same fibre with probability p = 1. Any two nodes which do not lie on the same fibre are connected with probability pðx; r; cÞ ¼ 1 where x is the distance between the nodes, r = 7 is an arbitrary steepness parameter, and c is a characteristic distance. Connecting nodes according to a distance dependent attachment function is a standard technique when studying percolation on spatial networks [12]. For better performance, only nodes separated by x < 2 are considered for connection. By varying the spatial coupling parameter c, the number of transverse connections between fibres can be controlled from a state where fibres are strongly connected at c � 1, to a state where fibres become increasingly disconnected as c ! 0. We can think of the coupling c as representing the grade of fibrosis. If the coupling is large, there is no, or only low-grade, interstitial fibrosis. As coupling is reduced, this is equivalent to a greater burden of interstitial fibrosis. Ideally, our distribution of spatial coupling should be modulated using individual-specific fibrosis data, however, this information is not available with the current datasets. For each dataset, spatial networks are generated with five values of the characteristic distance, c. The specific values used and corresponding risk parameters (see section 3.5.1), are given in S2 Table in SM section 1.3.1 in S1 File. For convenience, specific parameter choices are labelled as having a low, medium or high global risk of micro-reentry.

Constructing spatial networks without fibre structure
For comparison purposes, we consider a null model where the spatial network excludes fibre orientation. All other geometric information is retained. The fibre-less spatial network is generated analogously to the fibre model, with the omission of the initial fibre tractography steps, illustrated schematically in Fig 4. Nodes are randomly distributed across each atrial mesh at a density of approximately 1 node per voxel and are added to the spatial network if the distance between the new node and any existing nodes is greater than or equal to d sep = 0.7. To construct the spatial network, nodes are connected using the same distance dependent probability function as the fibre model given in Eq (1). Unlike the fibre model, none of the nodes are already connected by edges which represent the underlying fibre structure of the atria (no pink edges); the probability of connecting two nodes is independent of their orientation.
For each atrial dataset, we generate spatial networks for three values of the characteristic coupling, corresponding to high, medium and low risk cases, see S3 Table in SM section 1.3.2 in S1 File. The characteristic coupling values for the fibre-less null model are significantly larger than for the fibre model since there are no permanent longitudinal edges.

Identifying the substrate for micro-anatomical reentry
For both the fibre and fibre-less models, we associate the substrate for micro-anatomical reentry with regions in the corresponding spatial network where closed loops exist of length ℓ � τ.
Here, τ represents a typical refractory wavelength. These are found by applying a discrete diffusion model (DDM) to the network, see SM section 1.4.  (1), where c is a characteristic distance (dashed circle). There is no directional bias in the connection probability (blue edges only). Nodes which are closer together than the characteristic distance (green arrows) are more likely to be connected than those which are far apart (red arrows). (C) The micro-reentrant substrate is identified using the DDM. https://doi.org/10.1371/journal.pone.0267166.g004 The DDM is related to several extremely simple physics models of micro-anatomical reentry [19][20][21][22][23]. These models are not electrophysiologically realistic models of AF. The strength of discrete models lies in the ease with which structural discontinuities can be modelled, in contrast to conventional models which are significantly complicated when no-flux boundary conditions are imposed at the local level [24]. This justifies the use of the DDM in the current context where we focus exclusively on the structural basis for micro-reentry. However, this approach is not suitable for studying the dynamics of micro-anatomical reentry. For a wider discussion of the role of model choice see [25].

Defining risk of micro-anatomical reentry.
To derive a measure for the risk of micro-reentry in each spatial network, we use the DDM and assume that all nodes in the network are equally susceptible to conduction block. If each structure has a similar number of nodes which, if blocked, may initiate a micro-anatomical reentry, then the probability of initiating any given circuit is approximately constant. Hence, for a spatial network with N potential substrates, the rate, λ, at which new micro-reentries are initiated in the DDM is directly related to the number of substrates N. This argument follows directly from similar arguments in [21]. A mathematically precise formulation of this argument is provided in SM section 1.5. The rate parameter, λ, is used as our measure of the global risk of micro-anatomical reentry across the atria as a whole. A uniform probability of conduction block is likely an unrealistic assumption. However, realistic conduction block behaviour is difficult to obtain in the DDM framework. These issues are discussed in more detail in the limitations section.
To compute the risk of micro-anatomical in a given atrial region, we define the local risk for a specific voxel, v, by calculating the number of circuits which form in that particular voxel, N v , when sampling a fixed number of identified circuits, N. The rescaled rate of reentry for the voxel v is then given byR where the denominator provides a normalisation. If a specific atrial region is referred to as "at risk", this signifies a non-zero value of the local risk,R v , in that region. Note that summing the local risk over all the voxels in the atrial geometry gives the global risk, l ¼ P vRv . Given that there are typically many more voxels in each atrial geometry than sampled circuits, we smooth the local risk using a 3d Gaussian convolution with a standard deviation of 5 voxels (* 1.5mm) when plotting onto the atrial surfaces.

Spatial clustering of reentrant circuits.
For each spatial network, we measure how spatially clustered the identified micro-reentrant substrate is. To do so, we sample 1000 voxels from the set of voxels in which micro-anatomical reentry is identified for each network and extract the coordinate of each sampled voxel. The number of independent spatial clusters is then derived by using the DBSCAN algorithm, a standard algorithm for spatial data clustering, see [26]. We set eps = 10 (� 3mm), which is the parameter controlling the distance between points for these to be considered in the same cluster. The number of points required for a cluster is set to 1. If only one cluster is identified by DBSCAN, this indicates that all 1000 sampled circuits in a given network fall within a single, small region of the atrial structure. Conversely, if 1000 clusters are identified, no two circuits are identified in the same location.

Atrial wall thickness & occupied voxel fraction.
Wall thickness is a common measurement used when analysing atrial geometry. One issue when calculating wall thickness is that it is not easily defined in the atrial bulk. On the surfaces, wall thickness can be defined as the distance through the atrial wall along the perpendicular surface vector (there are competing, but similar definitions). However, in the bulk many perpendicular surface vectors pass through the same voxel so that the thickness is not uniquely defined.
To resolve this issue, we introduce a novel measure analogous to wall thickness, the occupied voxel fraction (OVF), which quantifies the proximity of voxels to the atrial walls and captures differences in the local wall curvature and thickness gradients, see SM section 1.6. The measure is defined as the number of voxels within a radiusr ¼ 5 (approximately 1.5mm) which are inside the atrial structure, normalised by the total number of voxels within a sphere of radiusr. The OVF increases (decreases) if the atrial walls become thicker (thinner), and distinguishes between regions with convex, flat, or concave wall morphologies.

Longitudinal connection fraction.
For the spatial networks with fibre structure, we would like to associate the identified micro-reentrant substrate with the degree of local fibre alignment in a particular atrial region. Fibre alignment is not naturally defined in a spatial network. Therefore, we use the longitudinal connection fraction (LCF) as a proxy measure. This is defined as the average number of longitudinal connections (along the fibre direction), normalised by the total number of connections, either longitudinal or transverse. If the LCF in a voxel is large, then fibres in that region are well aligned. Conversely, if the LCF is small, the fibre structure is highly disordered in a local area. Note, however, that the LCF can increase both from an increase in the number of longitudinal connections, but also from a reduction in the number of transverse connections. Therefore, larger LCF does not necessarily imply more longitudinal connections. More detail is provided in SM section 1.7 where we validate that the LCF in each spatial network accurately reflects the expected degree of fibre misalignment measured in the underlying data.
It is interesting to note that we can compare the LCF to the degree of microscopic anisotropy, calculated using Eq (1), at the highest level of spatial coupling where micro-anatomical reentry is inducible in our spatial networks. This demonstrates that the onset of micro-reentry in our spatial networks occurs at approximately the same degree of structural anisotropy as found experimentally by Spach et al. [27], see SM section 1.8 for details. Figs 5 & 6 show the substrate identified as susceptible to micro-reentry in the human atria using the fibre and fibre-less models respectively, at low, medium and high global risk (progressively larger reentry rate, λ). Note that low, medium and high global risk refers to the overall risk of micro-reentry across the whole atrial tissue, and does not refer to whether specific regions identified are at low or high risk (which is captured by the local risk,R v ). The spatial distributions of local risk for the sheep atria are shown in SM section 2.1. How the differences in our datasets affect our results is discussed in SM section 2.2.

Spatial distribution of identified micro-reentrant substrate
The reader should be conscious that the local risk distributions shown do not represent the regions at risk for a single atrial tissue, but rather show the collation of the local risk regions identified across 1000 randomly generated spatial networks. If a region is highlighted in the left column, the low global risk case, this may be interpreted as an area which is susceptible to micro-reentry even if the fibrosis density is low (but non-zero). Regions identified in the middle and right columns are regions susceptible to micro-reentry at moderate and high fibrosis densities respectively. This illustrates a key result: different regions of the atria are susceptible to micro-reentry at different fibrosis densities. In particular, regions may have a characteristic fibrosis density range within which they are susceptible to micro-reentry; if the fibrosis density is too high or low these regions will not exhibit micro-reentry. At high global risk, the susceptible substrate is significantly more diffuse than at lower global risk. However, there remain multiple regions with very low or zero local risk, in particular the posterior right atrium (PRA), left atrial appendage (LAA) and inferior vena cava (IVC).
The dominant local risk substrate lies along the strip between two small holes in the superior right atrium. Secondary regions of high local risk are maintained in the PVs, with local risk migrating slightly further from the LSPV junction into the PLA. The opening of the the coronary sinus forms a new substrate which is not observed at higher coupling. The local risk observed previously at the IAS and PMs is largely absent. 4.1.2 Micro-reentrant substrate in spatial network without fibre structure. The local risk substrate identified by the fibre-less null model, see Fig 6, is qualitatively similar to the local risk regions identified in the fibre model. Key local risk regions including the superior right atrium and the PV sleeves are common to both models. Likewise, the PRA, LAA and IVC are not susceptible to micro-reentry in either model.
At low global risk, the dominant local risk substrate in the fibre-less null model is confined to two spatially isolated PM ridges, one of which exhibits local risk in the fibre model. Additional diffuse local risk is observed in the superior right atrium and in the right PV sleeves. At medium global risk, the low local risk substrate is retained, with additional local risk in all four PVs and diffuse risk across the superior right atrium. These local risk regions are consolidated at high global risk, with a small reduction in local risk observed along the spatially isolated PV ridges. The substrates observed for the fibre model along the coronary sinus (CS) opening and at the IAS are not observed in the fibre-less null model.

Clustering of micro-reentrant substrate
The observed micro-reentrant risk substrates for the fibre and fibre-less models suggest that micro-reentry is spatially confined for both models at low global risk, but that the substrate covers a wider area as the global risk of micro-reentry increases. Fig 7 shows that the substrate is more clustered for the fibre model than for the fibre-less null model. One possible interpretation of this result is that if all micro-reentrant circuits are confined to a small number of clusters, these will be easy to destroy or isolate via catheter ablation. However, if circuits are widely distributed, they will be more difficult, if not impossible, to destroy or isolate.  significantly below the average value for each dataset. As global risk increases (increasing λ), the mean OVF increases. This indicates that the bias to thin convex regions of the atria is reduced as edges are progressively removed in each spatial network. Framed in terms of fibrosis densities, this result illustrates that regions with low OVF are susceptible to micro-reentry at a lower fibrosis density than regions with larger OVF. These results are supported by statistical analysis in SM section 2.3. In the superior right atrium, the micro-reentrant substrate at high global risk is concentrated along a strip between two small structural holes (region 1) for the fibre model, Fig 9(A). This region exhibits moderate risk in the fibre-less null model, Fig 9(B), with additional risk along the spatially isolated PM (region 2), the junction of the CT and a second PM ridge (region 3), and around the opening of the TV (region 4). Region 2 exhibits low local risk in the fibre model, despite being a low OVF region, but is the dominant risk region for the fibre-less null model when the rate of micro-reentry is low or medium, see Fig 6(E) left and centre columns. This is because the strong longitudinal coupling in region 2 suppresses micro-reentry in the fibre model. This effect is absent in the null model where fibre structure is omitted. Region 2 is most susceptible to micro-reentry in the fibre-less null model at a low and moderate fibrosis densities, with the local risk reducing at higher densities. Regions 3 and 4 show very low or zero local risk in the fibre model.

Comparing fibre and fibre-less micro-reentrant substrates
Comparing the identified risk regions to the OVF values shown in Fig 9(C), we note that regions 1-4 all exhibit low OVF values (blue on the colour scale). Of the four regions identified, we observe that the LCF values, Fig 9(D), in region 1 are low (blue), indicating fibre misalignment, whereas strong longitudinal coupling is observed in regions 2-4. This suggests that the local risk of micro-reentry is suppressed in the fibre model in regions with high LCF, but can be enhanced in regions with low LCF. Recall that the OVF is a geometric property of the atrial structure which is unaffected by the choice of coupling parameter, c, used to initialise the spatial networks. Hence, the distributions in Fig 10(A) & 10(D) align along the OVF axis. However, the LCF is not independent of the coupling parameter c used to initialise each network, and, consequently, is not independent of the global risk, λ. Therefore, the distributions in Fig 10(A) & 10(D) do not align along the LCF axis since an increase in global risk is associated with a decrease in transverse coupling, c. This is equivalent to an increase in the LCF.  Fig 10(B), the voxel distribution is significantly skewed for both the OVF and the LCF with the weighted mean falling in the first OVF decile and the first LCF quartile. For the fibre-less null model at low global risk, Fig 10(C), the mean OVF value also falls in the first decile, but there is a small skew to LCF values larger than the median. This is an indirect effect and is explained by the strong longitudinal coupling along spatially isolated regions such as the PMs.
For both the fibre and fibre-less models the local risk density lies almost exclusively in the first OVF quartile, with negligible density at higher OVF values. In contrast, LCF values are clearly skewed in the fibre model, although some local risk density remains at larger LCF values. This implies that at low global risk, the OVF is the dominant factor in determining the local risk substrate, with fibre structure playing an important secondary role.
As the global risk of micro-reentry increases, we observe that the bias to low OVF values remains, but is slightly reduced, see  LCF) has a protective effect, reducing the local risk of micro-reentry, relative to regions with lower LCF. However, if we compare the case of low global risk in Fig 10(A)-10(C), to the case of high global risk, Fig 10(D)-10(F), it appears the the increase in global risk is associated with a shift in the voxel distributions from low LCF to higher LCF. This would would seem to imply, incorrectly, that longitudinal connections enhance the local risk of micro-reentry.
To understand this, note that there are two ways in which the LCF may increase: Firstly, if we fix the number longitudinal connections in the atrial geometry and progressively remove transverse connections, this will result in the average LCF increasing since the relative abundance of longitudinal connections increases. Conversely, if we fix the number of transverse connections and increase the number of longitudinal connections, this will result in an increase in the LCF. The former is the mechanism relevant when comparing the low global risk case to the high global risk case, and is not associated with an increase in the number of longitudinal connections. However, the latter is the mechanism relevant in a given atrial geometry where transverse coupling is effectively fixed, but regions with higher LCF have more longitudinal connections than regions with lower LCF. It is this second mechanism that is associated with the protective effect of longitudinal connections.

Discussion
In this paper, we introduce a novel computational approach for analysing how the structural substrate for micro-anatomical reentry may develop as the atria accumulate diffuse interstitial fibrosis. The technique combines spatial networks with fibre tractography, and assumes that interstitial fibrosis accumulation can be modelled as the decoupling of a spatial network, a concept closely related to spatial percolation [12]. Although our work focuses on micro-anatomical reentry, we believe our approach is of general interest, offering new methods for modelling atrial fibre structures, as well as metrics for quantifying longitudinal coupling and structural curvature.
Many studies investigate the electro-anatomical basis for micro-reentry as a potential mechanism of AF, although robust clinical evidence to support the mechanism is still lacking. Broadly speaking, such studies fall into three categories. In the first, the substrate can be identified directly, given sufficiently high resolution imaging, by observing the presence of a microanatomical driver. Examples include [7,9], identifying several key factors in the formation of the reentrant substrate. Of particular note are the orientation of muscle fibres, the structure, thickness, and thickness gradients of the atria, and the accumulation of fibrosis, particularly interstitial fibrosis [28].
Attempting to quantify these factors, the authors in [9] applied optical mapping to explanted human atria to identify how specific values of wall thickness and fibrosis burden correlate to the location of the known micro-reentrant substrate [9]. Driver regions were found to correlate well with areas of 20-30% wall thickness and 20-30% fibrotic burden, notably the junction of the RIPV and PLA, and between the CT and fibrotically insulated PMs. Fibre misalignment was also implicated, specifically between the CT and the PMs which also exhibit abrupt changes in the local wall thickness. However, these studies took a largely static view of the micro-reentrant substrate, not considering how the substrate may change under electrical or structural remodelling. Note that these studies have also been criticised for inducing drivers under physiologically unrealistic doses of pinacidil for shortening the action potential duration.
In the second category, the substrate for micro-reentry has been established in the atria, but is hidden due to the lack of a clearly observable driver. One recent study has addressed this problem by noting that the visible micro-reentrant substrate varies strongly with variable atrial refractoriness [29]. Hence, the authors were able to demonstrate that the hidden micro-reentrant substrate could be unmasked and stabilised by shortening atrial refractoriness with adenosine. Such changes, induced pharmacologically, are analogous to some of the electromechanical changes that may be expected over time from atrial remodelling [4]. However, a robust framework for predicting the future micro-reentrant substrate is still lacking.
No study is yet to address the third category which is to predict how the substrate for micro-reentry will develop in the future on a patient-specific basis, a problem directly related to the increasing need for arrhythmic risk stratification [30].
In this proof of concept-which we acknowledge is far from clinical applicability-we have attempted to take a step towards addressing this and ask how the substrate for micro-reentry may develop over time. Specifically, we ask how different parts of the atria are susceptible to micro-reentry at different characteristic fibrosis ranges; some regions of the atria, particularly those where the atrial walls are thin and there is significant fibre misalignment, may be susceptible to micro-reentry at a low fibrotic density, whereas other regions may require higher fibrosis densities before micro-reentry can be induced. Such an approach is currently only possible in-silico, primarily due to ongoing challenges with accurate in-vivo atrial imaging [31].
Our results implicate the role of wall thickness and the misalignment of fibres, suggesting that thin atrial walls and reduced longitudinal fibre coupling both enhance the probability that a region is susceptible to micro-reentry at lower fibrosis densities, supporting the findings in [9]. However, our study suggests that the dependence on these factors evolves as the density of interstitial fibrosis grows in the atria. In particular, the spatial spread of the micro-reentrant substrate increases dramatically with small reductions in the spatial network coupling, and indicates that the bias to thin atrial regions with complex fibre morphology is reduced as the micro-reentrant substrate becomes more spatially diverse. Many of the specific regions highlighted as local risk substrates for micro-reentry in [9,29], such as PM ridges and the superior right atrium, naturally emerge as local risk substrates using our method, most likely due to their position in thin atrial regions.
One explanation for the importance of wall thickness, fibrosis density, and local wall curvature is that driver regions emerge from percolation-like dynamics where the critical fibrosis density is dependent on the thickness of the structure studied [32,33]. This in turn may result in micro-anatomical circuits anchoring to the atrial surfaces in paroxysmal AF but distributing across the atrial wall in persistent AF [34]. Aside from micro-reentry, modelling fibrosis using percolation-style distributions is known to perform particularly well when simulating patterns of AF maintenance [35], and has been used to explain complex fractionated electrograms and reentries in fibrotic border zones [36,37].
In the wider AF literature, a number of computational studies consider the role of atrial structure on AF dynamics. In most cases, these studies investigate how structural factors affect electrical wavefront dynamics by solving the mono-or bidomain equation coupled with a suitable atrial cell model [33,[38][39][40][41][42][43]. However, even the most advanced patient-specific modelling methodologies suffer from the continued struggle to extract precise fibre orientation data and high resolution fibrosis profiles (particularly diffuse interstitial fibrosis) in-vivo [6].
Until these issues are resolved, computational modelling must continue to be used as a platform for hypothesis testing, studying the dynamics of AF initiation and maintenance with a variety of methods and across scales. In particular, the clearest insights may be found from models which contrast the role of a specific electro-anatomical feature with null models in the absence of that feature. Examples include models with and without realistic tissue geometries [38], with and without patient-specific fibrosis [44], isotropic vs. anisotropic fibre structures [41], and continuous vs. discrete modelling choices [25]. We hope that the techniques introduced in this paper can add to the range of techniques employed for hypothesis testing in cardiac electrophysiology.

Micro-anatomical reentry in the wider context
The role of local drivers, and specifically micro-anatomical reentrant circuits, in initiating and maintaining AF is controversial. After initial findings showed promise in 2015 [7], some suggested that the mechanism may have "the potential to unify some of the previous discrepant observations" [8] on AF mechanisms. However, although research is ongoing, reproducible clinical evidence identifying micro-anatomical reentrant circuits in patients is still lacking. Despite this, it is important to reiterate that the size of micro-reentrant circuits are postulated to be at, or below, the spatial resolution which can be resolved with conventional multi-electrode mapping [7,9]. This may hamper attempts to acquire clinical evidence in support of the mechanism.
Opposing the local drivers hypothesis, many believe that AF is maintained by spatio-temporal chaos in the atria, whereby wavefronts continually collide into each other producing new fibrillatory wavefronts [45][46][47][48][49][50][51]; evidence for the mechanism is extensive and extends beyond AF-specific studies to more general studies on spatio-temporal chaos in excitable systems [52][53][54]. In our view, there is no reason that local drivers cannot co-exist with other non-local mechanisms, especially given recent evidence for a spectrum of mechanisms at different levels of organisation during fibrillation [4]. This is important because the interaction between local drivers and spatio-temporal chaos may have a significant effect on whether the ablation of local drivers is a feasible strategy for terminating AF.
Spatio-temporal chaos may have a non-trivial impact on local driver regions, likely resulting in an increased rate of circuit termination as the fibrotic substrate becomes more spread out, but also resulting in AF becoming more turbulent and difficult to terminate. If local drivers can be ablated, this may be avoided. However, local ablation in this context with numerous drivers across a diffuse substrate may become practically impossible given that (1) circuits are short-lived due to chaotic propagation, (2) circuits are spatially diverse, and (3) even if circuits can be destroyed, damage to the atrial myocardium may promote turbulent electrical propagation, driving more persistent AF. Despite this, local ablation strategies may still be worth pursuing if they lower the probability of initiating AF from sinus rhythm, where micro-reentrant circuits may act as a trigger.
Whether robust evidence for the micro-reentrant driver mechanism can be established is a question for the future. However, even without this evidence, the methods introduced here may be valuable to the wider cardiac electrophysiology community, adding to the computational toolkit for hypothesis testing. In particular, the methods introduced for constructing spatial networks using fibre orientation data may be of general interest to the community. Equally, structural metrics including the occupied voxel fraction and the longitudinal connection fraction may be useful for studying atrial geometries and fibre structures. Finally, although we acknowledge the limitations of the method (see below), the discrete diffusion model used to propagate signals across our spatial networks may be of general use, particularly in cases of high computational complexity where more phenomenologically accurate models constrain the number of simulations that can be run.

Limitations
The aims of this study are highly focused, discussing the structural basis for micro-anatomical reentry only. Many other factors may affect the probability that a micro-reentrant circuit forms in a given region of the atria, including other forms of fibrosis and ionic remodelling.
Our approach does not consider these factors and cannot exclude their importance to microanatomical reentry.
The key technical limitations of our work come under two main categories: (1) Imaging related limitations and (2) limitations arising from specific modelling choices in the construction of the spatial networks.
In our view, one of the most important limitations relates to the datasets used for analysis. In particular, the human fibre dataset has undergone extensive smoothing, lacking realistic local heterogeneity in the fibre structure; the observed heterogeneous regions are likely artefacts from the synthetic fibre generation method. Similarly, the atrial geometry is heavily smoothed, lacking fine structural detail in regions such as the RAA and LAA, but with numerous structural holes across the geometry. Any future work must ensure that predictions are not simply an artefact of the data acquisition and generation process.
Although using high resolution fibre maps like those in [16,17] may avoid these issues, acquiring such data is not feasible in-vivo, remaining a major ongoing research challenge [55,56].
Assuming the imaging data used is of an acceptable quality and accuracy, a number of limitations arise in the fibre map modelling process. Firstly, the fibre tractography methods applied to generate global fibre tracts from local fibre orientation data are only an approximation of reality. Without performing precise histology to determine the position of global fibre tracts in our model, not histology focusing solely on local fibre orientation, it is not possible for us to explicitly validate this approach. Since the imaging data used was acquired from previous studies, such a process is not possible in the current work but may be possible in future studies. As a simple check of the robustness of our results, the full simulation pipeline, including the regeneration of fibres from a different set of random seed points, was carried out three times for the healthy sheep atria, showing no evidence that this significantly affected our final results.
Once fibre tracts are generated, nodes are placed along the fibres, and are coupled to nearby nodes based on their separation. The attachment function to connect two nodes is uniform across the atrial structure and does not consider any possible differences in the local connectivity of cardiomyocytes. Likewise, the characteristic coupling, c, that is used in each model is constant across the atrial structure. This can be thought of as applying a uniform density of interstitial fibrosis across the atria. It is known that not all regions of the atria are equally susceptible to the accumulation of fibrosis [57], and likewise, different forms of fibrosis accumulate differently across the atria [58]. Our assumption of uniformity is a reasonable first approximation in the absence of patient-specific information regarding the distribution of interstitial fibrosis which may not be visible with LGE-MRI. We do not consider any macroscopic fibrosis in the model. Such information was not available with the datasets acquired, but should be considered in any future work. However, how to accurately map imagingderived fibrosis densities to precise spatial network edge densities is so far unclear and will likely pose a future technical challenge.
Once a spatial network has been constructed, the micro-reentrant substrate is identified by applying a simple discrete diffusion model on the network. The approach is loosely based on the techniques discussed in [20,34], although we stress that its purpose here is strictly to identify locations in which isolated fibres result in micro-anatomical reentry, not to simulate AF dynamics. Assessing the ability of these structures to initiate and maintain AF would warrant a more detailed computational analysis with a phenomenologically accurate model which may be challenging at the resolution of our spatial network, and may struggle to produce large scale statistics. Despite this, the use of phenomenologically accurate propagation models may be necessary within our spatial network framework if we are to fully gauge and weigh up the relative impact of structural factors on the probability of micro-reentry, relative to other considerations such as the formation of functional blocks. Future studies using either method would benefit from a detailed comparison to experimentally acquired atrial functional data.
The regions that are identified by the discrete diffusion model are a set of isolated fibre tracts where unidirectional conduction block at a single node is sufficient to induce a microreentrant circuit. The method to apply conduction block assumes uniform risk of block, avoiding the need to specify special electrical properties in key atrial locations, and ignoring the potential impact of functional blocks from source-sink mismatch. In practice, it is known that electrical properties of cardiomyocytes vary significantly across the atria. For instance, the proarrhythmic conditions found in the PVs such as conduction velocity slowing and shortened action potential durations are not considered. However, our results do indicate that the PVs can emerge as a key risk substrate without the inclusion of electrical proarrhythmic effects.
Finally, before our approach can have real clinical relevance, it would benefit from further validation. We have demonstrated that the spatial networks accurately preserve the underlying properties of each atrial fibre map, and that our techniques ensure an even density of nodes across the network. For clinical relevance, it is important to validate the approach directly against experimental data, ideally with datasets with a history of micro-anatomical reentry. If possible, raw data should be available such that fibre maps can be generated using a range of methods at different levels of smoothing, and that the influence of small structural holes can be tested. Data may be acquired at different levels of interstitial fibrotic density to assess the validity of our method's longitudinal predictions.

Conclusion
We have introduced a simple, proof of concept framework which attempts to study how the substrate for micro-anatomical reentry develops from the accumulation of interstitial fibrosis. The method, which is based on the application of percolation to spatial networks, suggests that the micro-reentrant substrate is critically dependent of local tissue geometry and areas of fibre misalignment. We suggest that the dependence on these factors is complex, continuously evolving with the absolute level of global micro-reentrant risk.
If robust clinical evidence can be found to support the micro-reentrant driver mechanism, and if our methods develop sufficiently, the approach introduced here may have potential for patient-specific risk stratification and personalised ablation strategies. However, our results also imply that the micro-reentrant substrate may become so spatially diverse, particularly at high fibrosis densities, that ablating these drivers may not be practically feasible.
Supporting information S1 File. Supplementary information document. This file contains additional details regarding the methods used in this paper, as well as expanded results. In particular, this document includes details of atrial geometry thickness measurements, and figures showing the microanatomical risk substrates for the sheep atrial datasets. (PDF)